In [4]:
Im
Out[4]:
Al final no hemos hecho esto
We import the math libraries, data analysis libraries and image libraries.
# math libraries
import numpy as np
import sympy
from sympy import *
t = Symbol('t')
from scipy.interpolate import interp2d
from scipy.interpolate import griddata
import math
# data analysis library
import pandas as pd
# import plotting libraries
import plotly.offline as go_offline
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.ticker as ticker
# imaging library
from PIL import Image
We set the data and imagen directories
DATADIR='data/' # Directory with the data
FIGURESDIR='figures/' # Figures produced
To obtain the topography of the study area we delimited the area and then used the Copernicus DEM tiff image eu_dem_v11_E30N10.Tiff. Since this image weighs 1.1 G, we crop it to get a smaller image that looks like
Im=Image.open(FIGURESDIR+'eu_dem_v11_E30N10_crop.png')
Im
We build a mesh, each point of this mesh is given by its UTM_X, UTM_Y and elevation. The data for this mesh is saved as a cvs file that we open.
# Topografy data
topo=pd.read_csv(DATADIR+'topografy.csv')
topo
| UTM_X | UTM_Y | elevation | |
|---|---|---|---|
| 0 | 618131.414840 | 4.197675e+06 | 702.294556 |
| 1 | 618132.694676 | 4.197586e+06 | 707.380493 |
| 2 | 618133.974489 | 4.197497e+06 | 722.496460 |
| 3 | 618135.254279 | 4.197409e+06 | 721.459595 |
| 4 | 618136.534045 | 4.197320e+06 | 732.357910 |
| ... | ... | ... | ... |
| 2585 | 623304.834524 | 4.195087e+06 | 1092.907593 |
| 2586 | 623306.169278 | 4.194999e+06 | 1152.172607 |
| 2587 | 623307.504009 | 4.194910e+06 | 1132.570068 |
| 2588 | 623308.838716 | 4.194821e+06 | 1062.827759 |
| 2589 | 623310.173398 | 4.194732e+06 | 986.483154 |
2590 rows × 3 columns
We split the data in topo to get three numpy arrays
topo1=topo[['UTM_X','UTM_Y','elevation']].to_numpy()
topox=[x[0] for x in topo1]
topoy=[x[1] for x in topo1]
topoelv=[x[2] for x in topo1]
In the list cotasxy we save the dimensions of the model.
cotasxy=[min(topox)+200,max(topox)-2000,min(topoy)+100,max(topoy)-200]
The following function surface interpolates the data on a grid of points, the interpolation method can be linear, cubic or nearest.
def surface(grid,metodo): #metodo='linear', 'cubic' o 'nearest'
tx=grid[0]
ty=grid[1]
tz=grid[2]
X = np.linspace(min(tx), max(tx))
Y = np.linspace(min(ty), max(ty))
X, Y = np.meshgrid(X, Y) # 2D grid for interpolation
Z = griddata((tx,ty),tz,(X,Y),
method=metodo)
return [tx,ty,tz,X,Y,Z]
We choose the linear interpolation method and apply it to the topographic data.
topo_gr=[topox,topoy,topoelv]
topo_linear=surface(topo_gr,'linear')
We are ready to create the figure that realizes the 3D model. We use two color scales 'YlGn' and 'gray' to get more depth in the model.
fig=go.Figure()
# topografía
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='YlGn',
opacity = 0.8,
name='Topography',
legendgroup='topo',
showlegend=True,
showscale=False))
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='gray',
opacity = 0.3,
name='Topo',
legendgroup='topo',
showlegend=False,
showscale=False))
camera = dict(up=dict(x=0, y=0, z=1),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=-1, z=2) )
fig.update_layout( title="Topography Model",
#paper_bgcolor = 'black',
scene = dict(
xaxis=dict(title='UTM_X',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[0],cotasxy[1]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
yaxis=dict(title='UTM_Y',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[2],cotasxy[3]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
zaxis=dict(nticks=4,
tickfont = dict(size = 10,color = 'black'),
title='Elevation',
title_font_size=10,
titlefont_color='black',
range=[500,1000],
backgroundcolor='white',
color='black',
gridcolor='gray'),
aspectratio=dict(x=2, y=1.8, z=0.3)),
#aspectmode='data'),
scene_camera= camera
)
fig.show()
With the following function contact_tr we process the data from the cvs files in which the stratigraphic and mechanical data are stored
def contact_tr(csv,aumento_elevacion):
tr=pd.read_csv(csv)
tr1=tr[['x','y','SAMPLE_1']].to_numpy()
tx=[x[0] for x in tr1]
ty=[x[1] for x in tr1]
telv=[x[2]+aumento_elevacion for x in tr1]
return [tx,ty,telv]
Load and process the data
# mechanic contacts ????
palomeque=contact_tr(DATADIR+'palomeque.csv',0)+['palomeque_thrust']
palomeque_ds=contact_tr(DATADIR+'palomeque_ds.csv',0)+['palomeque_desplazado']
calvillo=contact_tr(DATADIR+'calvillo.csv',0)+['calvillo_thrust']
calvillo_ds=contact_tr(DATADIR+'calvillo_ds.csv',0)+['calvillo_desplazado']
# faults
fault2=contact_tr(DATADIR+'fault2.csv',0)+['fault2']
fault1=contact_tr(DATADIR+'fault1.csv',0)+['fault1']
fault3=contact_tr(DATADIR+'fault3.csv',0)+['fault1']
# stratigraphic contacts ?????
E12cal=contact_tr(DATADIR+'E1E2 calvillo.csv',0)+['E1-E2 calvillo']
E12pal=contact_tr(DATADIR+'E1E2 palomeque.csv',0)+['E1-E2 palomeque']
E23pal=contact_tr(DATADIR+'E2E3 palomeque.csv',0)+['E2-E3 palomeque']
E3Opal=contact_tr(DATADIR+'E3O palomeque.csv',0)+['E3-O palomeque']
PEpal=contact_tr(DATADIR+'PE palomeque.csv',0)+['P-E palomeque']
PEcal=contact_tr(DATADIR+'PE calvillo.csv',0)+['P-E calvillo']
E12cal_ds=contact_tr(DATADIR+'E1E2 calvillo_ds.csv',0)+['E1-E2 calvillo ds']
E12pal_ds=contact_tr(DATADIR+'E1E2 palomeque_ds.csv',0)+['E1-E2 palomeque ds']
The following function contact_dat uses the comad Scatter3d to create de contact data for the figure.
def contact_dat(x,y,elv,h,mode,nam,color,gr,t,ls):
z=[x+h for x in elv]
return go.Scatter3d(x=x, y=y, z=z,
mode =mode,
name=nam,
legendgroup=gr,
showlegend=t,
line=dict(color=color,
dash=ls,
width=5)
)
cont=[palomeque_ds,palomeque,calvillo,calvillo_ds]
cab=[E12cal,E12pal,E23pal,E3Opal,PEpal,PEcal,E12cal_ds,E12pal_ds]
fal=[fault1,fault2,fault3]
contacts=[contact_dat(palomeque_ds[0],palomeque_ds[1],palomeque_ds[2],30,
"lines",'Thrusts','black','contacts',True,None)
]+[contact_dat(x[0],x[1],x[2],30,"lines",x[3],'black','contacts',False,None) for x in cont[1:]]
cabs=[contact_dat(E12cal[0],E12cal[1],E12cal[2],30,"lines",'Stratigraphic contacts','black',
'cabalgamientos',True,'dot')
]+[contact_dat(x[0],x[1],x[2],30,"lines",x[3],'black','cabalgamientos',False,'dot') for x in cab[1:]]
faults=[contact_dat(fault1[0],fault1[1],fault1[2],50,"lines",'Strike-slip faults','black','faults',True,None)
]+[contact_dat(x[0],x[1],x[2],30,"lines",x[3],'black','faults',False,None) for x in fal[1:]]
fig=go.Figure(contacts+cabs+faults)
# topography
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='YlGn',
opacity = 0.8,
name='Topography',
legendgroup='topo',
showlegend=True,
showscale=False))
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='gray',
opacity = 0.3,
name='Topo',
legendgroup='topo',
showlegend=False,
showscale=False))
camera = dict(up=dict(x=0, y=0, z=1),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=-1, z=2) )
fig.update_layout( title="Topography Model",
#paper_bgcolor = 'black',
scene = dict(
xaxis=dict(title='UTM_X',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[0],cotasxy[1]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
yaxis=dict(title='UTM_Y',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[2],cotasxy[3]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
zaxis=dict(nticks=4,
tickfont = dict(size = 10,color = 'black'),
title='Elevation',
title_font_size=10,
titlefont_color='black',
range=[500,1000],
backgroundcolor='white',
color='black',
gridcolor='gray'),
aspectratio=dict(x=2, y=1.8, z=0.3)),
#aspectmode='data'),
scene_camera= camera
)
fig.show()
We add texts for heights and Stratigraphic levels to the figure
namesx=[619200,620800]
namesy=[4196500,4196250]
namesz=[850,850]
namestxt=['Calvillo height','Palomeque height','A','A$^\prime$']
tx=[calvillo[0][1]+50,calvillo[0][3],calvillo[0][4],calvillo[0][5],calvillo[0][7]]
ty=[calvillo[1][1]-50,calvillo[1][3]-50,calvillo[1][4]-70,calvillo[1][5]-50,calvillo[1][7]-70]
tz=[calvillo[2][1]+40,calvillo[2][3]+30,calvillo[2][4]+50,calvillo[2][5]+15,calvillo[2][7]+20]
nnx=[618900,619350,620200,619350,619400,620120,620700,
618850,619200,620900,621100,619700]
nny=[4196800,4196700,4196400,4196200,4196000,4196100,4196000,
4195500,4195500,4195800,4195500,4195200]
nnz=[700,750,750,750,700,790,700,750,700,600,650,700]
nntxt=['M1','P','P','E1','E2','E1','E2','E1','E2','E3','O2','E1']
fig=go.Figure(contacts+cabs+faults)
# topography
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='YlGn',
opacity = 0.8,
name='Topography',
legendgroup='topo',
showlegend=True,
showscale=False))
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='gray',
opacity = 0.3,
name='Topo',
legendgroup='topo',
showlegend=False,
showscale=False))
trace = go.Scatter3d(
x=namesx, y=namesy, z=namesz,
text=namestxt,
name='Heights',
mode="text",
textfont=dict(color=["black","black"],size=13),
hoverinfo="skip")
fig.add_trace(trace)
trace2 = go.Scatter3d(
x=nnx, y=nny, z=nnz,
text=nntxt,
name='Stratigraphic levels',
mode="text",
legendgroup='cabalgamientos',
showlegend=False,
textfont=dict(color='black',size=18),
hoverinfo="skip")
fig.add_trace(trace2)
camera = dict(up=dict(x=0, y=0, z=1),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=-1, z=2) )
fig.update_layout( title="Topography Model",
#paper_bgcolor = 'black',
scene = dict(
xaxis=dict(title='UTM_X',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[0],cotasxy[1]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
yaxis=dict(title='UTM_Y',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[2],cotasxy[3]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
zaxis=dict(nticks=4,
tickfont = dict(size = 10,color = 'black'),
title='Elevation',
title_font_size=10,
titlefont_color='black',
range=[500,1000],
backgroundcolor='white',
color='black',
gridcolor='gray'),
aspectratio=dict(x=2, y=1.8, z=0.3)),
#aspectmode='data'),
scene_camera= camera
)
fig.show()
To add triangles to the contact tracks we need some mathematical functions
def eq(A,B): # A and B are points
x1, y1 = A
x2, y2 = B
# calculate the distance between the two points
distance = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
# calculate the angle between the two points
angle = math.atan2(y2 - y1, x2 - x1)
# calculate the coordinates of the third point (the first one)
x3_1 = x1 + distance * math.cos(angle - (1 * math.pi / 3))
y3_1 = y1 + distance * math.sin(angle - (1 * math.pi / 3))
# calculate the coordinates of the third point (the second one)
x3_2 = x1 + distance * math.cos(angle + (1 * math.pi / 3))
y3_2 = y1 + distance * math.sin(angle + (1 * math.pi / 3))
return[[x3_1,y3_1],[x3_2,y3_2]]
def tr(i,contacto,h,s,hh):
A=(contacto[0][i],contacto[1][i])
c1z=contacto[2][i]+h
c2x=contacto[0][i+1]
c2y=contacto[1][i+1]
d=(c2x-A[0],c2y-A[1])
nd=np.sqrt(d[0]**2+d[1]**2)
dn=(d[0]/nd,d[1]/nd)
B=(A[0]+dn[0]*150,A[1]+dn[1]*150)
C=eq(A, B)[s]
return [[A[0],B[0],C[0]],[A[1],B[1],C[1]],[c1z,c1z,c1z+hh]]
tcalvillo=[tr(1,calvillo,15,0,100),tr(3,calvillo,15,0,50),tr(5,calvillo,15,0,50),tr(7,calvillo,15,0,0)]
tcalvillo_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
alphahull=5, opacity=1, color='black',
i = np.array([0]),
j = np.array([1]),
k = np.array([2]),
legendgroup='contacts',
showlegend=False,
) for x in tcalvillo]
tcalvillo_ds=[tr(1,calvillo_ds,25,1,30)]
tcalvillo_ds_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
alphahull=5, opacity=1, color='black',
i = np.array([0]),
j = np.array([1]),
k = np.array([2]),
legendgroup='contacts'
) for x in tcalvillo_ds]
tpalomeque=[tr(1,palomeque,15,0,70),tr(4,palomeque,15,0,70),tr(5,palomeque,15,0,90),
tr(6,palomeque,15,0,70),tr(7,palomeque,15,0,30)]
tpalomeque_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
alphahull=5, opacity=1, color='black',
i = np.array([0]),
j = np.array([1]),
k = np.array([2]),
legendgroup='contacts'
) for x in tpalomeque]
tpalomeque_ds=[tr(1,palomeque_ds,5,1,0)]
tpalomeque_ds_dat=[go.Mesh3d(x=x[0],y=x[1],z=x[2],
alphahull=5, opacity=1, color='black',
i = np.array([0]),
j = np.array([1]),
k = np.array([2]),
legendgroup='contacts',
) for x in tpalomeque_ds]
The triangles data are stored in tcont.
tcont=tcalvillo_dat+tcalvillo_ds_dat+tpalomeque_dat+tpalomeque_ds_dat
fig=go.Figure(contacts+cabs+faults+tcont)
# topography
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='YlGn',
opacity = 0.8,
name='Topography',
legendgroup='topo',
showlegend=True,
showscale=False))
fig.add_trace(go.Surface( x=topo_linear[3],y=topo_linear[4],z=topo_linear[5],
colorscale='gray',
opacity = 0.3,
name='Topo',
legendgroup='topo',
showlegend=False,
showscale=False))
trace = go.Scatter3d(
x=namesx, y=namesy, z=namesz,
text=namestxt,
name='Heights',
mode="text",
textfont=dict(color=["black","black"],size=13),
hoverinfo="skip")
fig.add_trace(trace)
trace2 = go.Scatter3d(
x=nnx, y=nny, z=nnz,
text=nntxt,
name='Stratigraphic levels',
mode="text",
legendgroup='cabalgamientos',
showlegend=False,
textfont=dict(color='black',size=18),
hoverinfo="skip")
fig.add_trace(trace2)
camera = dict(up=dict(x=0, y=0, z=1),
center=dict(x=0, y=0, z=0),
eye=dict(x=0, y=-1, z=2) )
fig.update_layout( title="Topography Model",
#paper_bgcolor = 'black',
scene = dict(
xaxis=dict(title='UTM_X',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[0],cotasxy[1]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
yaxis=dict(title='UTM_Y',
tickfont = dict(size = 10,color = 'black'),
title_font_size=10,
titlefont_color='black',
range=[cotasxy[2],cotasxy[3]],
backgroundcolor='white',
color='black',
gridcolor='gray'),
zaxis=dict(nticks=4,
tickfont = dict(size = 10,color = 'black'),
title='Elevation',
title_font_size=10,
titlefont_color='black',
range=[500,1000],
backgroundcolor='white',
color='black',
gridcolor='gray'),
aspectratio=dict(x=2, y=1.8, z=0.3)),
#aspectmode='data'),
scene_camera= camera
)
fig.show()